/* search.c is part of VIP program
 * Written by J. Salvador Arias, 2010
 *      CONICET, INSUE, Instituto Miguel Lillo,
 *      Miguel Lillo 205, S.M. de Tucuman (4000)
 *      Tucuman, Argentina
 *      e-mail: jsalarias@csnat.unt.edu.ar
 *
 * VIP is a free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public
 * License as published by the Free Software Foundation and
 * Creative Commons.
 *
 * This file and the VIP program is distributed in the hope
 * that it will be useful, but WITHOUT ANY WARRANTY; without
 * even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE.
 * See the GNU General Public License for more details.
 *
 * To view a copy of this license, visit
 * http://creativecommons.org/licenses/GPL/2.0/
 * or send a letter to
 * Creative Commons, 171 Second Street, Suite 300,
 * San Francisco, California, 94105, USA.
 */

#include "biogeo.h"
#include "bitfield.h"
#include "grid.h"
#include "tessellation.h"

#define JUST_COUNT_COST     1
#define FULL_OPTIMIZATION   2
#define SEARCH_OPTIMIZATION 3

/* Code specificly used for searches and optimizations */

static SREAL OptimizeNode (GNode* node, SFLAG mode);
static void OptimizeAll (GPtrArray* recArray, SFLAG optimize);

static void CopyRec (RECONS* dest, RECONS* src);
static GPtrArray* CreateRecArrayCopy (GPtrArray* src);
static GPtrArray* CreateRecArrayCopyOnNode (GPtrArray* trgt, GPtrArray* src, GNode* node);
static void CopyRecArray (GPtrArray* dest, GPtrArray* src);

static SFLAG AreRecsEqual (RECONS* probe, RECONS* test);
static SFLAG CompareRec (GPtrArray* probe, GList* recList);
static SREAL CheckReconstruction (GList* recs, GPtrArray* toCheck);
static SFLAG AddReconstruction (GList** recs, GPtrArray* toAdd, SFLAG check, SFLAG copy);

static SREAL IncrementalDownPass (GNode* startNode);

static void MeasureOverlap (RECONS* actRec);

static SREAL NodeRemoval (PHYLOTREE* actTree, GPtrArray* srcRec, GPtrArray* sectorNodeArray, SCOUNTER hold, GList** recs);
static SREAL BlindAnnealing (PHYLOTREE* actTree, GPtrArray* recArray, PROGRESSSTRUCT* progress);
static SREAL GibbsAnnealing (PHYLOTREE* actTree, GPtrArray* recArray, GPtrArray* sectorArray, INTARRAY* prob, PROGRESSSTRUCT* progress);
static SREAL NixonAnnealing (PHYLOTREE* actTree, GPtrArray* recArray, GPtrArray* sectorArray, PROGRESSSTRUCT* progress);

static void UpdateRecCounts (GPtrArray* recArray);

static void SetDefaultBremer (GPtrArray* recArray);
static void UpdateBremer (GNode* startNode, SREAL target);
static void SetBremer (GPtrArray* recArray);

static SFLAG IsRemovalSecure (PHYLONODE* toCheck);

static GPtrArray* tmpNodeArray = NULL;

void SetRec (GPtrArray* recArray) {
 SCOUNTER i;
 RECONS* actRec;

 for (i = recArray -> len - 1; i >= 0; -- i) {
   actRec = (RECONS*) g_ptr_array_index (recArray, i);
   actRec -> node -> rec = actRec;
 }
}

void MakeOptimization (SFLAG justCount) {
 GList* navTree;
 GList* navRec;
 PHYLOTREE* tree;
 GNode* navNode;
 PHYLONODE* actNode;
 RECONS* actRec;
 GPtrArray* rec;

 for (navTree = trees; navTree != NULL; navTree = g_list_next (navTree)) {
   tree = (PHYLOTREE*) navTree -> data;
   if (tree -> recOR == NULL) {
     /* There are no reconstruction, first creates the array */
     tree -> recOR = g_ptr_array_new ();
     for (navNode = tree -> root; navNode != NULL; ) {
       actNode = (PHYLONODE*) navNode -> data;
       UpdateOnCalculus ();
       if (actNode -> rec == NULL) {
         actRec = CreateRecons (actNode);
         g_ptr_array_add (tree -> recOR, actRec);
         actNode -> rec = actRec;
         if (g_node_first_child (navNode) != NULL) {
           navNode = g_node_first_child (navNode);
           continue;
         }
       }
       if (g_node_next_sibling (navNode) == NULL)
         navNode = navNode -> parent;
       else
         navNode = g_node_next_sibling (navNode);
     }
     actNode = (PHYLONODE*) tree -> root -> data;
     OptimizeAll (tree -> recOR, TRUE);
     tree -> best = tree -> worst = actNode -> rec -> cost;
     SetDefaultBremer (tree -> recOR);
   }
   else {
     /* If there are a reconstruction, updates the costs of the reconstructions */
     SetRec (tree -> recOR);
     if (justCount == TRUE)
       OptimizeAll (tree -> recOR, FALSE);
     else
       OptimizeAll (tree -> recOR, TRUE);
     actNode = (PHYLONODE*) tree -> root -> data;
     tree -> best = tree -> worst = actNode -> rec -> cost;
     SetDefaultBremer (tree -> recOR);
     if (tree -> recs != NULL) {
       for (navRec = tree -> recs; navRec != NULL; navRec = g_list_next (navRec)) {
         UpdateOnCalculus ();
         rec = (GPtrArray*) navRec -> data;
         SetRec (rec);
         if (justCount == TRUE)
           OptimizeAll (rec, FALSE);
         else
           OptimizeAll (rec, TRUE);
         SetBremer (rec);
         actNode = (PHYLONODE*) tree -> root -> data;
         if (actNode -> rec -> cost < tree -> best)
           tree -> best = actNode -> rec -> cost;
         if (actNode -> rec -> cost > tree -> worst)
           tree -> worst = tree -> worst;
       }
     }
     SetRec (tree -> recOR);
   }
 }
}

SREAL OptimizeNode (GNode* node, SFLAG mode) {
 RECONS* actRec;
 PHYLONODE* actNode;
 PHYLONODE* descNode;
 PHYLONODE* left;
 PHYLONODE* right;
 GNode* navDesc;
 SREAL timeFraq;
 SINTEGER ageDif;

 actNode = (PHYLONODE*) node -> data;
 actRec = actNode -> rec;
 actRec -> cost = 0;
 actRec -> numVics = actRec -> numNoVics = actRec -> numRemovNodes = 0;
 actRec -> nodeFlag &= REMOV_NODE;        /* If the node is a removed, then nodeFlag is equal to REMOV_NODE, otherwise nodeFlag is 0 */
 if (actRec -> triangles != NULL)
   DestroyTriangleList (actRec -> triangles);
 actRec -> triangles = NULL;
 if (g_node_first_child (actNode -> node) == NULL) {
   if (mode == FULL_OPTIMIZATION) {
     CopyBitField (actRec -> obs, actNode -> term -> obs);
     if (CountBits (actRec -> obs) > 0) {
       DoFill (actRec -> cells, actRec -> obs, actRec -> node -> fill);
       actRec -> numCells = CountBits (actRec -> cells);
       actRec -> numObs = CountBits (actRec -> obs);
       actRec -> maxAge = actNode -> minAge;
     }
     else {         /* An empty set "optimization" */
       ResetBitField (actRec -> cells);
       actRec -> numCells = actRec -> numObs = 0;
     }
   }
   if (actRec -> numCells == 0) return 0;   /* Empty set exit */
   if (actRec -> nodeFlag & REMOV_NODE) {
     actRec -> numRemovNodes = 1;
     actRec -> cost = params . costRemov;
   }
   return actRec -> cost;
 }
 if (mode & (FULL_OPTIMIZATION | SEARCH_OPTIMIZATION)) {
   ResetBitField (actRec -> obs);
   ResetBitField (actRec -> cells);
   actRec -> numObs = 0;
   actRec -> numCells = 0;
   actRec -> active = 0;
   actRec -> maxAge = actNode -> minAge;
   actRec -> setA = NULL;
   actRec -> setB = NULL;
   if (!(params . search & USE_LATEOVER))
     actRec -> removSet = NULL;
 }
 left = right = NULL;
 for (navDesc = g_node_first_child (actNode -> node); navDesc != NULL; navDesc = g_node_next_sibling (navDesc)) {
   descNode = (PHYLONODE*) navDesc -> data;
   actRec -> numVics += descNode -> rec -> numVics;
   actRec -> numNoVics += descNode -> rec -> numNoVics;
   actRec -> numRemovNodes += descNode -> rec -> numRemovNodes;
   actRec -> cost += descNode -> rec -> cost;
   if (mode & (FULL_OPTIMIZATION | SEARCH_OPTIMIZATION)) {
     if ((descNode -> rec -> numCells == 0) || (descNode -> rec -> nodeFlag & REMOV_NODE)) continue;
     if (descNode != actRec -> removSet) {
       BitFieldORsum (actRec -> obs, descNode -> rec -> obs);
       BitFieldORsum (actRec -> cells, descNode -> rec -> cells);
       if ((actRec -> maxAge != NO_AGE) && (actRec -> maxAge > descNode -> rec -> maxAge) && (descNode -> rec -> maxAge != NO_AGE))
         actRec -> maxAge = descNode -> rec -> maxAge;
     }
     ++ actRec -> active;
     if (left == NULL)
       left = descNode;
     else if (right == NULL)
       right = descNode;
   }
 }
 if (mode & (FULL_OPTIMIZATION | SEARCH_OPTIMIZATION)) {
   actRec -> numCells = CountBits (actRec -> cells);
   if (actRec -> numCells == 0) {
     /* This is an empty set */
     actRec -> numObs = 0;
     actRec -> cost = 0;
     return actRec -> cost;
   }
   actRec -> numObs = CountBits (actRec -> obs);
 }
 if (actRec -> nodeFlag & REMOV_NODE) {
   ++ actRec -> numRemovNodes;
   actRec -> cost += params . costRemov;
 }
 if (actRec -> removSet != NULL) {
   ++ actRec -> numRemovNodes;
   actRec -> cost += params . costPartRemov;
 }
 if (actRec -> active == 1) {
   actRec -> nodeFlag &= ~ VIC_NODE;
   if (params . costRemov == 0) {
     ++ actRec -> numNoVics;
     actRec -> cost += params . costNoVic;
   }
   if ((params . search & CONSTR_SEARCH) && (actNode -> constrained))
     actRec -> cost += 10000;
   return actRec -> cost;
 }
 else if (actRec -> active != 2) {
   actRec -> nodeFlag &= ~ VIC_NODE;
   ++ actRec -> numNoVics;
   actRec -> cost += params . costNoVic;
   if ((params . search & CONSTR_SEARCH) && (actNode -> constrained))
     actRec -> cost += 10000;
   return actRec -> cost;
 }
 if (mode & (FULL_OPTIMIZATION | SEARCH_OPTIMIZATION)) {
   actRec -> timeLag = actNode -> minAge - actRec -> maxAge;
   /* Preserves the order, indep of tree sorting */
   if (left -> id < right -> id) {
     actRec -> setA = left;
     actRec -> setB = right;
   }
   else {
     actRec -> setB = left;
     actRec -> setA = right;
   }
   MeasureOverlap (actRec);
 }
 if ((actRec -> overlap <= params . overlap) && (actRec -> overlap >= 0)) {
   if (params . mode & USE_TIME) {
     ageDif = actRec -> timeLag - params . timeLag;
     if (ageDif > 0) {
       if ((params . mode & TIME_REAL) && (ageDif < params . timeDiff)) {
         timeFraq = ((SREAL) ageDif / params . timeDiff) * params . costTimeLag;
         actRec -> cost += timeFraq;
       }
       else
         actRec -> cost += params . costTimeLag;
       actRec -> nodeFlag |= TIME_LAG;
     }
   }
   /* Reset the age of the node */
   actRec -> maxAge = actNode -> minAge;
   actRec -> nodeFlag |= VIC_NODE;
   ++ actRec -> numVics;
   if (params . mode & OVER_REAL)
     actRec -> cost += actRec -> overlap;
 }
 else {
   actRec -> nodeFlag &= ~ VIC_NODE;
   ++ actRec -> numNoVics;
   actRec -> cost += params . costNoVic;
   if ((params . search & CONSTR_SEARCH) && (actNode -> constrained))
     actRec -> cost += 10000;
 }
 return actRec -> cost;
}

void OptimizeAll (GPtrArray* recArray, SFLAG optimize) {
 SINTEGER i;
 RECONS* actRec;

 for (i = recArray -> len - 1; i >= 0; -- i) {
   actRec = (RECONS*) g_ptr_array_index (recArray, i);
   UpdateOnCalculus ();
   OptimizeNode (actRec -> node -> node, FULL_OPTIMIZATION);
 }
}

void SetDefaultBremer (GPtrArray* recArray) {
 SREAL orVal;
 SINTEGER i;
 RECONS* actRec;

 actRec = (RECONS*) g_ptr_array_index (recArray, 0);
 orVal = actRec -> cost;
 for (i = 0; i < recArray -> len; ++ i) {
   actRec = (RECONS*) g_ptr_array_index (recArray, i);
   if (actRec -> nodeFlag & VIC_NODE)
     actRec -> node -> bremer = G_MAXDOUBLE;
   else
     actRec -> node -> bremer = orVal;
 }
}

void SetBremer (GPtrArray* recArray) {
 SREAL actVal;
 SINTEGER i;
 RECONS* actRec;

 actRec = (RECONS*) g_ptr_array_index (recArray, 0);
 actVal = actRec -> cost;
 for (i = 0; i < recArray -> len; ++ i) {
   actRec = (RECONS*) g_ptr_array_index (recArray, i);
   if (actRec -> nodeFlag & VIC_NODE) continue;
   if (actRec -> node -> bremer > actVal)
     actRec -> node -> bremer = actVal;
 }
}

void UpdateBremer (GNode* startNode, SREAL target) {
 GNode* navNode;
 PHYLONODE* actNode;

 for (navNode = startNode; navNode != NULL; navNode = navNode -> parent) {
   actNode = (PHYLONODE*) navNode -> data;
   if (actNode -> rec -> nodeFlag & VIC_NODE) continue;
   if (actNode -> bremer > target)
     actNode -> bremer = target;
 }
 if (actNode -> rec -> cost < actNode -> tree -> best)
   actNode -> tree -> best = actNode -> rec -> cost;
 if (actNode -> rec -> cost > actNode -> tree -> worst)
   actNode -> tree -> worst = actNode -> tree -> worst;
}

SREAL IncrementalDownPass (GNode* startNode) {
 GNode* navNode;
 SREAL toRetCost;

 toRetCost = G_MAXDOUBLE;
 for (navNode = startNode; navNode != NULL; navNode = navNode -> parent)
   toRetCost = OptimizeNode (navNode, SEARCH_OPTIMIZATION);
 UpdateBremer (startNode, toRetCost);
 return toRetCost;
}

void MeasureOverlap (RECONS* actRec) {
 SREAL divA, divB, maxDiv;
 GNode* navDesc;
 PHYLONODE* descNode;

 actRec -> overlap = 0;
 divA = (SREAL) BitFieldCommon (actRec -> setB -> rec -> cells, actRec -> setA -> rec -> obs) / actRec -> setA -> rec -> numObs;
 divB = (SREAL) BitFieldCommon (actRec -> setA -> rec -> cells, actRec -> setB -> rec -> obs) / actRec -> setB -> rec -> numObs;
 actRec -> overlap = (divA > divB) ? divA : divB;
 /* Check if there is no overlap between a node and the descendants of its sister
    This is to remove the effect of some spurious disjunction, which produce a double detection of
    the same barrier. As is overlap against a more deepest node, it "accepts" up to the double of
    the params . overlap setting */
 if ((actRec -> overlap > 0) && (actRec -> overlap <= params . overlap)) {
   for (navDesc = g_node_first_child (actRec -> setA -> node); navDesc != NULL; navDesc = g_node_next_sibling (navDesc)) {
     descNode = (PHYLONODE*) navDesc -> data;
     if ((descNode -> rec -> numCells == 0) || (descNode -> rec -> nodeFlag & REMOV_NODE)) continue;
     if (descNode == actRec -> setA -> rec -> removSet) continue;
     divA = (SREAL) BitFieldCommon (actRec -> setB -> rec -> cells, descNode -> rec -> obs) / descNode -> rec -> numObs;
     divB = (SREAL) BitFieldCommon (descNode -> rec -> cells, actRec -> setB -> rec -> obs) / actRec -> setB -> rec -> numObs;
     maxDiv = ((divA > divB) ? divA : divB) / 2;
     if (maxDiv > actRec -> overlap)
       actRec -> overlap = maxDiv;
   }
   for (navDesc = g_node_first_child (actRec -> setB -> node); navDesc != NULL; navDesc = g_node_next_sibling (navDesc)) {
     descNode = (PHYLONODE*) navDesc -> data;
     if ((descNode -> rec -> numCells == 0) || (descNode -> rec -> nodeFlag & REMOV_NODE)) continue;
     if (descNode == actRec -> setB -> rec -> removSet) continue;
     divB = (SREAL) BitFieldCommon (actRec -> setA -> rec -> cells, descNode -> rec -> obs) / descNode -> rec -> numObs;
     divA = (SREAL) BitFieldCommon (descNode -> rec -> cells, actRec -> setA -> rec -> obs) / actRec -> setA -> rec -> numObs;
     maxDiv = ((divA > divB) ? divA : divB) / 2;
     if (maxDiv > actRec -> overlap)
       actRec -> overlap = maxDiv;
   }
 }
}

SREAL NodeRemoval (PHYLOTREE* actTree, GPtrArray* srcRec, GPtrArray* sectorNodeArray, SCOUNTER hold, GList** recs) {
 SREAL bestScore, actScore;
 SFLAG foundBetter, foundEqual, foundOverlap;
 SFLAG removNode, partRemov;
 SFLAG tmpNodeFlag;
 SFLAG update, shuffle;
 GNode* navNode;
 PHYLONODE* actNode;
 PHYLONODE* tmpNode;
 PHYLONODE* tmpRemov;
 RECONS* actRec;
 SCOUNTER k, bestPos, nodeLen;
 SCOUNTER randSel;
 SCOUNTER numHold;

 foundBetter = TRUE;
 numHold = 0;
 if (recs == NULL)
   numHold = hold;
 nodeLen = sectorNodeArray -> len;
 shuffle = TRUE;
 while (foundBetter == TRUE) {
   foundBetter = FALSE;
   foundOverlap = FALSE;
   foundEqual = FALSE;
   if (shuffle == TRUE)
     ShufflePtrArray (sectorNodeArray);
   actNode = (PHYLONODE*) actTree -> root -> data;
   bestScore = actNode -> rec -> cost;
   k = bestPos = 0;
   update = TRUE;
   removNode = partRemov = FALSE;
   for (k = 0; k < sectorNodeArray -> len; ++ k) {
     navNode = (GNode*) g_ptr_array_index (sectorNodeArray, k);
     if (navNode -> parent == NULL) continue;
     actNode = (PHYLONODE*) navNode -> data;
     actRec = actNode -> rec;
     if (actRec -> numCells == 0) continue;
     tmpNode = (PHYLONODE*) navNode -> parent -> data;
     tmpRemov = NULL;
     if (((tmpNode -> rec -> removSet == actNode) || (actRec -> nodeFlag & REMOV_NODE)) &&
         (params . algorithm & FLIP_NODES)) {
       /* Make the backup of the node */
       tmpNodeFlag = actRec -> nodeFlag;
       /* Make the change */
       if (actRec -> nodeFlag & REMOV_NODE)
         actRec -> nodeFlag &= ~REMOV_NODE;
       else {
         tmpRemov = tmpNode -> rec -> removSet;
         tmpNode -> rec -> removSet = NULL;
       }
       actScore = IncrementalDownPass (navNode);
       if (actScore < bestScore) {
         foundBetter = TRUE;
         foundEqual = FALSE;
         bestScore = actScore;
         bestPos = k;
         removNode = partRemov = FALSE;
         if (params . algorithm & FAST_SEARCH) {
           update = FALSE;
           break;
         }
       }
       else if ((actScore <= (bestScore + params . extra)) && (foundBetter != TRUE) && (foundEqual != TRUE)) {
         bestPos = k;
         removNode = partRemov = FALSE;
         foundEqual = TRUE;
       }
     }
     else {
       if (IsRemovalSecure (tmpNode) == FALSE) continue;
       /* Make the backup of the node */
       tmpNodeFlag = actRec -> nodeFlag;

       randSel = 0;
       if (params . search & USE_LATEOVER)
         randSel = g_random_int_range (0, 100);
       if ((params . search & USE_LATEOVER) && (randSel >= 50)) {
         tmpNode -> rec -> removSet = actNode;
         actScore = IncrementalDownPass (navNode);
         if (actScore < bestScore) {
           foundBetter = TRUE;
           foundEqual = FALSE;
           bestScore = actScore;
           bestPos = k;
           partRemov = TRUE;
           removNode = FALSE;
           if (params . algorithm & FAST_SEARCH) {
             update = FALSE;
             break;
           }
         }
         else if ((actScore <= (bestScore + params . extra)) && (foundBetter != TRUE) && (foundEqual != TRUE)) {
           bestPos = k;
           partRemov = TRUE;
           removNode = FALSE;
           foundEqual = TRUE;
         }
         tmpNode -> rec -> removSet = NULL;
       }
       actRec -> nodeFlag |= REMOV_NODE;
       actScore = IncrementalDownPass (navNode);
       if (actScore < bestScore) {
         foundBetter = TRUE;
         foundEqual = FALSE;
         bestScore = actScore;
         bestPos = k;
         removNode = TRUE;
         partRemov = FALSE;
         if (params . algorithm & FAST_SEARCH) {
           update = FALSE;
           break;
         }
       }
       else if ((actScore <= (bestScore + params . extra)) && (foundBetter != TRUE) && (foundEqual != TRUE)) {
         bestPos = k;
         removNode = TRUE;
         partRemov = FALSE;
         foundEqual = TRUE;
       }
       if ((params . search & USE_LATEOVER) && (randSel < 50)) {
         tmpNode -> rec -> removSet = actNode;
         actRec -> nodeFlag = tmpNodeFlag;
         actScore = IncrementalDownPass (navNode);
         if (actScore < bestScore) {
           foundBetter = TRUE;
           foundEqual = FALSE;
           bestScore = actScore;
           bestPos = k;
           partRemov = TRUE;
           removNode = FALSE;
           if (params . algorithm & FAST_SEARCH) {
             update = FALSE;
             break;
           }
         }
         else if ((actScore <= (bestScore + params . extra)) && (foundBetter != TRUE) && (foundEqual != TRUE)) {
           bestPos = k;
           partRemov = TRUE;
           removNode = FALSE;
           foundEqual = TRUE;
         }
         tmpNode -> rec -> removSet = NULL;
       }
     }
     /* Restores the reconstruction */
     actRec -> nodeFlag = tmpNodeFlag;
     if (tmpRemov != NULL)
       tmpNode -> rec -> removSet = tmpRemov;
     IncrementalDownPass (navNode);
   }
   if ((foundEqual == TRUE) && (numHold < hold))
     foundBetter = TRUE;
   if (foundBetter == TRUE) {
     if (update == TRUE) {
       navNode = g_ptr_array_index (sectorNodeArray, bestPos);
       actNode = (PHYLONODE*) navNode -> data;
       actRec = actNode -> rec;
       if (removNode == TRUE)
         actRec -> nodeFlag |= REMOV_NODE;
       else if (partRemov == TRUE) {
         tmpNode = (PHYLONODE*) navNode -> parent -> data;
         tmpNode -> rec -> removSet = actNode;
       }
       else {
         if (actRec -> nodeFlag & REMOV_NODE)
           actRec -> nodeFlag &= ~ REMOV_NODE;
         else {
           tmpNode = (PHYLONODE*) navNode -> parent -> data;
           tmpNode -> rec -> removSet = NULL;
         }
       }
       bestScore = IncrementalDownPass (navNode);
     }
     if (foundEqual == TRUE) {
       g_ptr_array_remove_index_fast (sectorNodeArray, bestPos);
       shuffle = FALSE;
     }
     else {
       numHold = 0;
       if (recs == NULL)
         numHold = hold;
       shuffle = TRUE;
       sectorNodeArray -> len = nodeLen;
     }
     if ((recs != NULL) && (AddReconstruction (recs, srcRec, TRUE, TRUE) == TRUE))
       ++ numHold;
   }
 }
 sectorNodeArray -> len = nodeLen;
 return bestScore;
}

SREAL BlindAnnealing (PHYLOTREE* actTree, GPtrArray* recArray, PROGRESSSTRUCT* progress) {
 SREAL actScore, actCost, prevScore, diff, bestCost;
 SCOUNTER actGen, toSel, object;
 GNode* navNode;
 PHYLONODE* rootNode;
 PHYLONODE* removSet;
 SFLAG tmpNodeFlag;
 RECONS* actRec;
 GPtrArray* nuRec;

 rootNode = (PHYLONODE*) actTree -> root -> data;
 bestCost = 2 * rootNode -> rec -> cost;
 for (actGen = 0; actGen < params . generations; ++ actGen) {
   UpdateOnCalculus ();
   if (progress -> isActive == FALSE) break;
   prevScore = rootNode -> rec -> cost;
   toSel = g_random_int_range (0, tmpNodeArray -> len);
   navNode = (GNode*) g_ptr_array_index (tmpNodeArray, toSel);
   if (navNode -> parent == NULL) continue;
   /* Make the backup of the node */
   tmpNodeFlag = actRec -> nodeFlag;
   removSet = actRec -> removSet;
   /* Make the change */
   if (ChangeNode (navNode) == FALSE) continue;
   actScore = rootNode -> rec -> cost;
   diff = actScore - prevScore;
   if (diff >= 0)
     object = (params . acceptEqual * 10) / (diff + 1);
   else {
     if (CheckReconstruction (actTree -> recs, recArray) <= 0) {
       nuRec = CreateRecArrayCopy (recArray);
       SetRec (nuRec);
       actCost = NodeRemoval (actTree, nuRec, tmpNodeArray, 0, NULL);
       SetRec (recArray);
       if (AddReconstruction (&actTree -> recs, nuRec, TRUE, FALSE) == TRUE) {
         if (bestCost > actCost)
           bestCost = actCost;
       }
       else
         DestroySingleRec (nuRec);
     }
     object = 1100;
   }
   toSel = g_random_int_range (0, 1000);
   if (toSel >= object) {
     /* Restores the reconstruction */
     actRec -> nodeFlag = tmpNodeFlag;
     actRec -> removSet = removSet;
     IncrementalDownPass (navNode);
   }
 }
 return bestCost;
}

SREAL GibbsAnnealing (PHYLOTREE* actTree, GPtrArray* recArray, GPtrArray* sectorArray, INTARRAY* prob, PROGRESSSTRUCT* progress) {
 SREAL actScore, prevScore, diff, actCost, bestCost;
 SCOUNTER actGen, toSel, object, sumProb, i;
 GNode* navNode;
 PHYLONODE* actNode;
 PHYLONODE* rootNode;
 PHYLONODE* removSet;
 SFLAG tmpNodeFlag;
 RECONS* actRec;
 GPtrArray* nuRec;

 rootNode = (PHYLONODE*) actTree -> root -> data;
 bestCost = 2 * rootNode -> rec -> cost;
 for (actGen = 0; actGen < params . generations; ++ actGen) {
   UpdateOnCalculus ();
   if (progress -> isActive == FALSE) break;
   prevScore = rootNode -> rec -> cost;
   prob -> len = 0;
   sumProb = 0;
   /* Make the node list */
   if (sectorArray -> len > 0)
     g_ptr_array_remove_range (sectorArray, 0, sectorArray -> len);
   for (i = 0; i < params . sectSize; ++ i) {
     toSel = g_random_int_range (0, tmpNodeArray -> len);
     navNode = (GNode*) g_ptr_array_index (tmpNodeArray, toSel);
     if (navNode -> parent == NULL) continue;
     actRec = actNode -> rec;
     /* Make the backup of the node */
     tmpNodeFlag = actRec -> nodeFlag;
     removSet = actRec -> removSet;
     /* Make the change */
     if (ChangeNode (navNode) == FALSE) continue;
     actScore = rootNode -> rec -> cost;
     diff = actScore - prevScore;
     if (diff >= 0)
       object = (params . acceptEqual * 10) / (diff + 1);
     else
       object = 1000 + ((params . acceptEqual * 10) / (diff - 1));
     AddToArray (prob, object);
     g_ptr_array_add (sectorArray, navNode);
     sumProb += object;
     /* Restores the reconstruction */
     actRec -> nodeFlag = tmpNodeFlag;
     actRec -> removSet = removSet;
     IncrementalDownPass (navNode);
   }
   toSel = g_random_int_range (0, sumProb);
   diff = 0;
   for (i = 0; i < prob -> len; ++ i) {
     diff += GetFromArray (prob, i);
     if (diff > toSel) break;
   }
   navNode = (GNode*) g_ptr_array_index (sectorArray, i);
   ChangeNode (navNode);
   if (CheckReconstruction (actTree -> recs, recArray) <= 0) {
     nuRec = CreateRecArrayCopy (recArray);
     SetRec (nuRec);
     actCost = NodeRemoval (actTree, nuRec, tmpNodeArray, 0, NULL);
     SetRec (recArray);
     if (AddReconstruction (&actTree -> recs, nuRec, TRUE, FALSE) == TRUE) {
       if (bestCost > actCost)
         bestCost = actCost;
     }
     else
       DestroySingleRec (nuRec);
   }
 }
 return bestCost;
}

SREAL NixonAnnealing (PHYLOTREE* actTree, GPtrArray* recArray, GPtrArray* sectorArray, PROGRESSSTRUCT* progress) {
 SCOUNTER actGen, toSel, i;
 GPtrArray* nuRec;
 SREAL bestCost, actCost;
 GNode* navNode;
 PHYLONODE* actNode;
 PHYLONODE* rootNode;

 RECONS* actRec;

 rootNode = (PHYLONODE*) actTree -> root -> data;
 bestCost = 2 * rootNode -> rec -> cost;
 for (actGen = 0; actGen < params . generations; ++ actGen) {
   UpdateOnCalculus ();
   if (progress -> isActive == FALSE) break;
   if (tmpNodeArray -> len > 0)
     g_ptr_array_remove_range (tmpNodeArray, 0, tmpNodeArray -> len);
   CopyPtrArray (tmpNodeArray, actTree -> nodeArray);
   /* Make the node list */
   if (sectorArray -> len > 0)
     g_ptr_array_remove_range (sectorArray, 0, sectorArray -> len);
   for (i = 0; i < params . sectSize; ++ i) {
     toSel = g_random_int_range (0, tmpNodeArray -> len);
     navNode = (GNode*) g_ptr_array_index (tmpNodeArray, toSel);
     if (navNode -> parent == NULL) continue;
     /* Check if elimination is valid */
     actNode = (PHYLONODE*) navNode -> data;
     actRec = actNode -> rec;
     if (IsRemovalSecure (actNode) == FALSE) continue;
     SetNodeAsRemoved (navNode);
     g_ptr_array_add (sectorArray, navNode);
     /* Remove element from node list */
     g_ptr_array_remove_index_fast (tmpNodeArray, toSel);
   }
   NodeRemoval (actTree, recArray, tmpNodeArray, 0, NULL);
   /* Restore the nodes */
   for (i = 0; i < sectorArray -> len; ++ i) {
     navNode = (GNode*) g_ptr_array_index (sectorArray, i);
     ChangeNode (navNode);
   }
   if (tmpNodeArray -> len > 0)
     g_ptr_array_remove_range (tmpNodeArray, 0, tmpNodeArray -> len);
   CopyPtrArray (tmpNodeArray, actTree -> nodeArray);
   actCost = NodeRemoval (actTree, recArray, tmpNodeArray, 0, NULL);
   if (CheckReconstruction (actTree -> recs, recArray) <= 0) {
     nuRec = CreateRecArrayCopy (recArray);
     if (AddReconstruction (&actTree -> recs, nuRec, TRUE, FALSE) == TRUE) {
       if (bestCost > actCost)
         bestCost = actCost;
     }
     else
       DestroySingleRec (nuRec);
   }
 }
 return bestCost;
}

SREAL FuseRecons (PHYLOTREE* actTree, PROGRESSSTRUCT* progress) {
 GList* navOrigArray;
 GList* navTestArray;
 GPtrArray* origArray;
 GPtrArray* testArray;
 SCOUNTER i;
 GNode* parent;
 RECONS* testRec;
 RECONS* origRec;
 GList* nuList;
 GPtrArray* nuRec;
 SREAL actCost, bestCost;

 nuList = NULL;
 bestCost = GetBestScore (actTree -> recs);
 if (actTree -> stored < 2) return bestCost;
 for (navOrigArray = actTree -> recs; (navOrigArray != NULL) && (g_list_next (navOrigArray) != NULL); navOrigArray = g_list_next (navOrigArray)) {
   origArray = (GPtrArray*) navOrigArray -> data;
   SetRec (origArray);
   for (navTestArray = g_list_next (navOrigArray); navTestArray != NULL; navTestArray = g_list_next (navTestArray)) {
     UpdateOnCalculus ();
     if (progress -> isActive == FALSE) break;
     testArray = (GPtrArray*) navTestArray -> data;
     for (i = 0; i < testArray -> len; ++ i) {
       UpdateOnCalculus ();
       if (progress -> isActive == FALSE) break;
       testRec = (RECONS*) g_ptr_array_index (testArray, i);
       if (testRec -> node -> node == actTree -> root) continue;
       origRec = (RECONS*) g_ptr_array_index (origArray, i);
       if (AreRecsEqual (testRec, origRec) == TRUE) continue;
       if ((testRec -> nodeFlag & REMOV_NODE) && (!(origRec -> nodeFlag & REMOV_NODE))) {
         if (IsRemovalSecure (origRec -> node) == FALSE) continue;
       }
       origRec -> node -> rec = testRec;
       parent = origRec -> node -> node -> parent;
       actCost = IncrementalDownPass (parent);
       if (actCost <= bestCost) {
         nuRec = CreateRecArrayCopyOnNode (origArray, testArray, testRec -> node -> node);
         SetRec (nuRec);
         actCost = NodeRemoval (actTree, nuRec, tmpNodeArray, 0, NULL);
         SetRec (origArray);
         if (AddReconstruction (&nuList, nuRec, TRUE, FALSE) != TRUE)
           DestroySingleRec (nuRec);
         if (bestCost > actCost)
           bestCost = actCost;
       }
       origRec -> node -> rec = origRec;
       IncrementalDownPass (parent);
     }
   }
 }
 /* There are new reconstructions! */
 if (nuList != NULL) {
   for (navTestArray = nuList; navTestArray != NULL; navTestArray = g_list_next (navTestArray)) {
     nuRec = (GPtrArray*) navTestArray -> data;
     if (AddReconstruction (&actTree -> recs, nuRec, TRUE, FALSE) != TRUE)
       DestroySingleRec (nuRec);
   }
   g_list_free (nuList);
 }
 return bestCost;
}

void CopyRec (RECONS* dest, RECONS* src) {
 dest -> node = src -> node;
 CopyBitField (dest -> cells, src -> cells);
 CopyBitField (dest -> obs, src -> obs);
 dest -> setA = src -> setA;
 dest -> setB = src -> setB;
 dest -> removSet = src -> removSet;
 dest -> nodeFlag = src -> nodeFlag;
 dest -> numVics = src -> numVics;
 dest -> numRemovNodes = src -> numRemovNodes;
 dest -> numNoVics = src -> numNoVics;
 dest -> cost = src -> cost;
 dest -> numCells = src -> numCells;
 dest -> numObs = src -> numObs;
 dest -> active = src -> active;
 dest -> overlap = src -> overlap;
 dest -> maxAge = src -> maxAge;
 dest -> timeLag = src -> timeLag;
 dest -> triangles = NULL;
}

GPtrArray* CreateRecArrayCopyOnNode (GPtrArray* orig, GPtrArray* toAdd, GNode* node) {
 SCOUNTER i;
 GPtrArray* dest;
 RECONS* nuRec;
 RECONS* actRec;

 dest = g_ptr_array_sized_new (orig -> len);
 for (i = 0; i < orig -> len; ++ i) {
   nuRec = (RECONS*) g_slice_alloc (sizeof (RECONS));
   nuRec -> cells = CreateBitField ();
   nuRec -> obs = CreateBitField ();
   nuRec -> triangles = NULL;
   actRec = (RECONS*) g_ptr_array_index (orig, i);
   if (IsDesc (actRec -> node -> node, node) == TRUE)
     CopyRec (nuRec, (RECONS*) g_ptr_array_index (toAdd, i));
   else
     CopyRec (nuRec, (RECONS*) g_ptr_array_index (orig, i));
   g_ptr_array_add (dest, nuRec);
 }
 return dest;
}

GPtrArray* CreateRecArrayCopy (GPtrArray* src) {
 SCOUNTER i;
 GPtrArray* dest;
 RECONS* nuRec;

 dest = g_ptr_array_sized_new (src -> len);
 for (i = 0; i < src -> len; ++ i) {
   nuRec = (RECONS*) g_slice_alloc (sizeof (RECONS));
   nuRec -> cells = CreateBitField ();
   nuRec -> obs = CreateBitField ();
   nuRec -> triangles = NULL;
   CopyRec (nuRec, (RECONS*) g_ptr_array_index (src, i));
   g_ptr_array_add (dest, nuRec);
 }
 return dest;
}

void CopyRecArray (GPtrArray* dest, GPtrArray* src) {
 SCOUNTER i;

 for (i = 0; i < src -> len; ++ i)
   CopyRec ((RECONS*) g_ptr_array_index (dest, i), (RECONS*) g_ptr_array_index (src, i));
}

SFLAG AreRecsEqual (RECONS* probe, RECONS* test) {
 if (probe -> node != test -> node) return FALSE;
 if (g_node_first_child (probe -> node -> node) == NULL) return (probe -> nodeFlag == test -> nodeFlag);
 if (probe -> cost != test -> cost) return FALSE;
 if (probe -> nodeFlag != test -> nodeFlag) return FALSE;
 if (probe -> numVics != test -> numVics) return FALSE;
 if (probe -> numRemovNodes != test -> numRemovNodes)  return FALSE;
 if (probe -> active != test -> active) return FALSE;
 if (probe -> nodeFlag & VIC_NODE) {
   if (probe -> setA != test -> setA) return FALSE;
   if (probe -> setB != test -> setB) return FALSE;
   if (probe -> removSet != test -> removSet) return FALSE;
   if (probe -> overlap != test -> overlap) return FALSE;
 }
 if (probe -> numCells != test -> numCells) return FALSE;
 if (probe -> numObs != test -> numObs) return FALSE;
 if (params . mode & USE_TIME) {
  if (probe -> timeLag != test -> timeLag) return FALSE;
  if (probe -> maxAge != test -> maxAge) return FALSE;
 }
 if (BitFieldEqual (probe -> obs, test -> obs) != TRUE) return FALSE;
 if (BitFieldEqual (probe -> cells, test -> cells) != TRUE) return FALSE;
 return TRUE;
}

SFLAG CompareRec (GPtrArray* probe, GList* recList) {
 GList* navList;
 GPtrArray* recArray;
 SFLAG isEqual;
 RECONS* probeRec;
 RECONS* actRec;
 SCOUNTER i;

 for (navList = recList; navList != NULL; navList = g_list_next (navList)) {
   recArray = (GPtrArray*) navList -> data;
   isEqual = TRUE;
   for (i = 0; i < probe -> len; ++ i) {
     probeRec = (RECONS*) g_ptr_array_index (probe, i);
     actRec = (RECONS*) g_ptr_array_index (recArray, i);
     if (AreRecsEqual (probeRec, actRec) != TRUE) {
       isEqual = FALSE;
       break;
     }
   }
   if (isEqual == TRUE) return TRUE;
 }
 return FALSE;
}

SREAL CheckReconstruction (GList* recs, GPtrArray* toCheck) {
 RECONS* actRec;
 SREAL best;

 if (recs == NULL) return -1;   /* The checked rec is "best" than no rec! */
 best = GetBestScore (recs);
 actRec = (RECONS*) g_ptr_array_index (toCheck, 0);
 return actRec -> cost - best;
}

SFLAG AddReconstruction (GList** recs, GPtrArray* toAdd, SFLAG check, SFLAG copy) {
 SREAL diff;
 SFLAG filter;
 GPtrArray* recToAdd;
 GList* lastLink;

 filter = FALSE;
 if ((check == TRUE) && (*recs != NULL)) {
   diff = CheckReconstruction (*recs, toAdd);
   if (diff < params . extra)
     filter = TRUE;         /* A better reconstruction is found! */
   else if (diff > params . extra) return FALSE;            /* Is worst, so, it is not added */
 }
 /* Check is the rec to store is different from stored recs */
 if ((*recs != NULL) && (filter != TRUE)) {
   if (CompareRec (toAdd, *recs) == TRUE) return FALSE;
 }
 if (copy == TRUE)
   recToAdd = CreateRecArrayCopy (toAdd);
 else
   recToAdd = toAdd;
 lastLink = InsertListElement (NULL, recToAdd);
 if (*recs == NULL)
   *recs = lastLink;
 else
   PutAtListEnd (*recs, lastLink);
 if (filter == TRUE)
   FilterRecList (recs, -1, params . extra);
 return TRUE;
}

void SearchMain (PROGRESSSTRUCT* progress) {
 GList* navTree;
 PHYLOTREE* tree;
 GPtrArray* sectorArray;
 GPtrArray* recArray;
 SCOUNTER stored, bestHits;
 SCOUNTER i, j, k;
 INTARRAY* prob;
 SREAL actScore, bestScore;
 SFLAG algorithm;

 progress -> actBar1 = j = 0;
 prob = NULL;
 if (params . algorithm & USE_SECTORS)
   sectorArray = g_ptr_array_sized_new (params . sectSize);
 else
   sectorArray = NULL;
 if ((params . algorithm & GIBBS_ANNEAL) || (params . algorithm & NIXON_ANNEAL)) {
   if (sectorArray == NULL)
     sectorArray = g_ptr_array_sized_new (params . sectSize);
   prob = CreateIntArray (params . sectSize);
 }
 for (navTree = GETTREES; navTree != NULL; navTree = g_list_next (navTree)) {
   if (progress -> isActive == FALSE) break;
   tree = (PHYLOTREE*) navTree -> data;
   tmpNodeArray = g_ptr_array_sized_new (tree -> nodeArray -> len);
   if (tree -> recs == NULL)
     AddReconstruction (&tree -> recs, tree -> recOR, TRUE, TRUE);
   bestScore = GetBestScore (tree -> recs);
   stored = g_list_length (tree -> recs);
   bestHits = 0;
   recArray = CreateRecArrayCopy (tree -> recOR);
   for (i = 0; i < params . its; ++ i) {
     if (progress -> isActive == FALSE) break;
     progress -> actBar2 = i;
     UpdateOnCalculus ();
     if (tmpNodeArray -> len > 0)
       g_ptr_array_remove_range (tmpNodeArray, 0, tmpNodeArray -> len);
     CopyPtrArray (tmpNodeArray, tree -> nodeArray);
     CopyRecArray (recArray, tree -> recOR);
     SetRec (recArray);
     if (params . algorithm & USE_SECTORS) {
       if (params . algorithm & FIRST_SECTOR) {
         algorithm = params . algorithm;
         params . algorithm &= ~PAGE_SEARCH;
         params . algorithm |= FAST_SEARCH;
       }
       for (k = 0; k < (2 * tmpNodeArray -> len); k += params . sectSize) {
         ShufflePtrArray (tmpNodeArray);
         if (sectorArray -> len > 0)
           g_ptr_array_remove_range (sectorArray, 0, sectorArray -> len);
         CopyPtrArrayPart (sectorArray, tmpNodeArray, 0, params . sectSize);
         actScore = NodeRemoval (tree, recArray, sectorArray, 0, NULL);
         if (params . algorithm & FIRST_SECTOR) {
           params . algorithm = algorithm;
           break;
         }
       }
     }
     if (params . algorithm & BLIND_ANNEAL) {
       params . algorithm &= ~PAGE_SEARCH;
       params . algorithm |= FAST_SEARCH;
       NodeRemoval (tree, recArray, tmpNodeArray, 0, NULL);
       params . algorithm = algorithm;
       actScore = BlindAnnealing (tree, recArray, progress);
     }
     else if (params . algorithm & GIBBS_ANNEAL) {
       params . algorithm &= ~PAGE_SEARCH;
       params . algorithm |= FAST_SEARCH;
       NodeRemoval (tree, recArray, tmpNodeArray, 0, NULL);
       params . algorithm = algorithm;
       actScore = GibbsAnnealing (tree, recArray, sectorArray, prob, progress);
     }
     else if (params . algorithm & NIXON_ANNEAL) {
       params . algorithm &= ~PAGE_SEARCH;
       params . algorithm |= FAST_SEARCH;
       NodeRemoval (tree, recArray, tmpNodeArray, 0, NULL);
       params . algorithm = algorithm;
       actScore = NixonAnnealing (tree, recArray, sectorArray, progress);
     }
     else {
       actScore = NodeRemoval (tree, recArray, tmpNodeArray, params . hold, &tree -> recs);
       if (params . keepAll == TRUE)
         AddReconstruction (&tree -> recs, recArray, FALSE, TRUE);
     }
     if (bestScore > actScore) {
       tree -> best = bestScore = actScore;
       bestHits = 1;
       stored = 0;
     }
     else if (bestScore == actScore)
       ++ bestHits;
   }
   if (progress -> isActive == FALSE)
     g_string_append_printf (outLog,"\n\n**** Search interrupted ****");
   stored = g_list_length (tree -> recs) - stored;
   g_string_append_printf (outLog,"\n\nTree %ld (%s)\nBest score: %0.4f\nHits: %ld of %ld iterations\nStored reconstructions: %ld",
                                   j, tree -> name, tree -> best, bestHits, i, stored);
   tree -> stored = g_list_length (tree -> recs);
   SetRec (tree -> recOR);
   tree -> actRec = NULL;
   DestroySingleRec (recArray);
   g_ptr_array_free (tmpNodeArray, TRUE);
   ++ j;
   progress -> actBar1 = j;
 }
 if (sectorArray != NULL)
   g_ptr_array_free (sectorArray, TRUE);
 if (prob != NULL)
   DestroyIntArray (prob);
}

void RefineRecs (PROGRESSSTRUCT* progress) {
 GList* navTree;
 PHYLOTREE* tree;
 GList* navRec;
 GPtrArray* recArray;
 SREAL bestScore, actScore;
 SCOUNTER stored, j;

 progress -> actBar1 = j = 0;
 for (navTree = GETTREES; navTree != NULL; navTree = g_list_next (navTree)) {
   if (progress -> isActive == FALSE) break;
   tree = (PHYLOTREE*) navTree -> data;
   if (tree -> recs == NULL) continue;
   tmpNodeArray = g_ptr_array_sized_new (tree -> nodeArray -> len);
   CopyPtrArray (tmpNodeArray, tree -> nodeArray);
   bestScore = GetBestScore (tree -> recs);
   stored = g_list_length (tree -> recs);
   if (params . algorithm & REC_FUSING) {
     actScore = FuseRecons (tree, progress);
     if (bestScore > actScore) {
       tree -> best = bestScore = actScore;
       stored = 0;
     }
   }
   else {
     recArray = CreateRecArrayCopy (tree -> recOR);
     for (navRec = tree -> recs; navRec != NULL;) {
       if (progress -> isActive == FALSE) break;
       UpdateOnCalculus ();
       CopyRecArray (recArray, (GPtrArray*)  navRec -> data);
       SetRec (recArray);
       actScore = NodeRemoval (tree, recArray, tmpNodeArray, params . maxHold, &tree -> recs);
       if (bestScore > actScore) {
         tree -> best = bestScore = actScore;
         stored = 0;
         navRec = tree -> recs;
         continue;
       }
       navRec = g_list_next (navRec);
     }
     DestroySingleRec (recArray);
   }
   if (progress -> isActive == FALSE)
     g_string_append_printf (outLog,"\n\n**** Search interrupted ****");
   stored = g_list_length (tree -> recs) - stored;
   g_string_append_printf (outLog,"\n\nTree %ld (%s)\nBest score: %0.4f\nStored reconstructions: %ld",
                                   j, tree -> name, tree -> best, stored);
   tree -> stored = g_list_length (tree -> recs);
   SetRec (tree -> recOR);
   tree -> actRec = NULL;
   g_ptr_array_free (tmpNodeArray, TRUE);
   ++ j;
   progress -> actBar1 = j;
 }
}

SFLAG IsRemovalSecure (PHYLONODE* toCheck) {
 RECONS* actRec;
 PHYLONODE* tmpNode;

 if (toCheck -> node == toCheck -> tree -> root) return FALSE;
 actRec = toCheck -> rec;
 if (actRec -> numCells == 0) return FALSE;
 if (actRec -> nodeFlag & REMOV_NODE) return FALSE;
 tmpNode = (PHYLONODE*) toCheck -> node -> parent -> data;
 if ((tmpNode -> rec -> nodeFlag & REMOV_NODE) || (tmpNode -> rec -> active < 2)) return FALSE;
 if (g_node_first_child (toCheck -> node) != NULL) {
   if (actRec -> active < 2) return FALSE;
 }
 if (tmpNode -> rec -> removSet != NULL) {
   if (tmpNode -> rec -> removSet == toCheck) return FALSE;
   if (tmpNode -> rec -> active == 2) return FALSE;
 }
 return TRUE;
}

void SetNodeAsRemoved (GNode* node) {
 PHYLONODE* actNode;
 RECONS* actRec;

 actNode = (PHYLONODE*) node -> data;
 actRec = actNode -> rec;
 actRec -> nodeFlag |= REMOV_NODE;
 IncrementalDownPass (node);
}

void SetNodeAsPartialRemov (GNode* node) {
 PHYLONODE* ancNode;

 if ((params . search & USE_LATEOVER) && (node -> parent != NULL)) {
   ancNode = (PHYLONODE*) node -> parent -> data;
   ancNode -> rec -> removSet = (PHYLONODE*) node -> data;
   IncrementalDownPass (node);
 }
}

SFLAG ChangeNode (GNode* node) {
 PHYLONODE* actNode;
 PHYLONODE* ancNode;
 RECONS* actRec;

 actNode = (PHYLONODE*) node -> data;
 if (actNode -> tree -> root == node) return FALSE;
 actRec = actNode -> rec;
 ancNode = (PHYLONODE*) node -> parent -> data;
 if (actRec -> nodeFlag & REMOV_NODE) {
   actRec -> nodeFlag &= ~ REMOV_NODE;
   IncrementalDownPass (node);
 }
 else if ((params . search & USE_LATEOVER) && (ancNode -> rec -> removSet == actNode)) {
   ancNode -> rec -> removSet = NULL;
   SetNodeAsRemoved (node);
 }
 else if (IsRemovalSecure (actNode) == TRUE) {
   if ((params . search & USE_LATEOVER) && (ancNode -> rec -> nodeFlag & VIC_NODE))
     SetNodeAsPartialRemov (node);
   else
     SetNodeAsRemoved (node);
 }
 else return FALSE;
 return TRUE;
}

void UpdateRecCounts (GPtrArray* recArray) {
 PHYLONODE* actNode;
 PHYLONODE* descNode;
 RECONS* actRec;
 GNode* navDesc;
 SCOUNTER i;

 for (i = recArray -> len - 1; i >= 0; -- i) {
   actRec = (RECONS*) g_ptr_array_index (recArray, i);
   actNode = actRec -> node;
   actRec -> cost = 0;
   actRec -> numVics = actRec -> numNoVics = 0;
   if (g_node_first_child (actNode -> node == NULL)) continue;
   for (navDesc = g_node_first_child (actNode -> node); navDesc != NULL; navDesc = g_node_next_sibling (navDesc)) {
     descNode = (PHYLONODE*) navDesc -> data;
     actRec -> numVics += descNode -> rec -> numVics;
     actRec -> numNoVics += descNode -> rec -> numNoVics;
   }
   if (actRec -> nodeFlag & VIC_NODE)
     ++ actRec -> numVics;
   else
     ++ actRec -> numNoVics;
 }
}

void MakeConsensus () {
 GList* navTree;
 GList* navRec;
 GPtrArray* tmpRec;
 SCOUNTER i;
 PHYLOTREE* tree;
 RECONS* actRec;
 RECONS* srcRec;

 for (navTree = GETTREES; navTree != NULL; navTree = g_list_next (navTree)) {
   tree = (PHYLOTREE*) navTree -> data;
   if (tree -> recCons != NULL) {
     SetRec (tree -> recCons);
     continue;
   }
   tmpRec = (GPtrArray*) tree -> recs -> data;
   tree -> recCons = CreateRecArrayCopy (tmpRec);
   for (navRec = g_list_next (tree -> recs); navRec != NULL; navRec = g_list_next (navRec)) {
     tmpRec = (GPtrArray*) navRec -> data;
     for (i = 0; i < tmpRec -> len; ++ i) {
       actRec = (RECONS*) g_ptr_array_index (tree -> recCons, i);
       srcRec = (RECONS*) g_ptr_array_index (tmpRec, i);
       BitFieldANDsum (actRec -> cells, srcRec -> cells);
       BitFieldANDsum (actRec -> obs, srcRec -> obs);
       if ((actRec -> setA != srcRec -> setA) || (actRec -> setB != srcRec -> setB)) {
         actRec -> setA = NULL;
         actRec -> setB = NULL;
         actRec -> nodeFlag &= ~ VIC_NODE;
       }
       if (actRec -> removSet != srcRec -> removSet)
         actRec -> removSet = NULL;
       if ((actRec -> nodeFlag & VIC_NODE) && (!(srcRec -> nodeFlag & VIC_NODE)))
         actRec -> nodeFlag &= ~ VIC_NODE;
       if ((actRec -> nodeFlag & REMOV_NODE) && (!(srcRec -> nodeFlag & REMOV_NODE)))
         actRec -> nodeFlag &= ~ REMOV_NODE;
       if (actRec -> numRemovNodes < srcRec -> numRemovNodes)
         actRec -> numRemovNodes = srcRec -> numRemovNodes;
     }
   }
   tree -> suppNodes = 0;
   SetRec (tree -> recCons);
   UpdateRecCounts (tree -> recCons);
 }
}

GList* CopyActRec (PHYLOTREE* tree) {
 GPtrArray* recArray;
 GList* toRet;

 if (tree -> actRec == NULL)
   recArray = CreateRecArrayCopy (tree -> recOR);
 else
   recArray = CreateRecArrayCopy ((GPtrArray*) tree -> actRec -> data);
 toRet = InsertListElement (NULL, recArray);
 if (tree -> recs == NULL)
   tree -> recs = toRet;
 else
   PutAtListEnd (tree -> recs, toRet);
 ++ tree -> stored;
 return toRet;
}

void SetConstrains (SFLAG toCons) {
 PHYLOTREE* tree;
 GList* navTree;
 RECONS* actRec;
 SCOUNTER i;

 for (navTree = GETTREES; navTree != NULL; navTree = g_list_next (navTree)) {
   tree = (PHYLOTREE*) navTree -> data;
   for (i = tree -> recOR -> len - 1; i >= 0; -- i) {
     actRec = (RECONS*) g_ptr_array_index (tree -> recOR, i);
     if ((toCons == TRUE) && (actRec -> nodeFlag & VIC_NODE))
       actRec -> node -> constrained = TRUE;
     else
       actRec -> node -> constrained = FALSE;
   }
 }
}
